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Abstract. Univariate polynomial root-finding is a classical subject, still 
important for modern computing. Frequently one seeks just the real roots 
of a real coefficient polynomial. They can be approximated at a low com¬ 
putational cost if the polynomial has no nonreal roots, but for high degree 
polynomials, nonreal roots are typically much more numerous than the 
real ones. The challenge is known for long time, and the subject has been 
intensively studied. The Boolean cost bounds for the refinement of the 
simple and isolated real roots have been decreased to nearly optimal, 
but the success has been more limited at the stage of the isolation of 
real roots. We obtain substantial progress by applying the algorithm of 
IS82I for the approximation of the root radii, that is, the distances of the 
roots to the origin. Namely we isolate the simple and well conditioned 
real roots of a polynomial at the Boolean cost dominated by the nearly 
optimal bounds for the refinement of such roots. We also extend our 
algorithm to the isolation of complex, possibly multiple, roots and root 
clusters staying within the same (nearly optimal) asymptotic Boolean 
cost bound. Our numerical tests with benchmark polynomials performed 
with the IEEE standard double precision show that our nearly optimal 
real root-finder is practically promising. Our techniques are simple, and 
their power and application range may increase in combination with the 
known efficient methods. 

Keywords: Polynomials; Root-finding; Real root-finding; Root isolation; Root 

radii 

1 Introduction 

Assume a univariate polynomial of degree n with real coefficients, 

n n 

p(x) ='^2pi xl = PnY]_{x - x j), Pn¥=®i 

i=0 i=1 


(1) 


which has r real roots x \,... ,x r and s = (n — r )/2 pairs of non-real complex 
conjugate roots. In many applications, e.g., to algebraic and geometric optimiza¬ 
tion, one seeks only the real roots, which make up just a small fraction of all 
roots. This motivates a well studied subject of real root-finding (see, e.g., [PT131 . 
[PT151 , [SM151 . and the extensive bibliography therein), but the most popular 
packages of subroutines for numerical root-finding such as MPSolve 2.0 IBF00] . 
Eigensolve IF021 . and MPSolve 3.0 [BR14] still approximate the r real roots 
about as fast and as slow as all the n complex roots. 

A typical fast real root-finder consists of two stages. At first one isolates all 
simple and well conditioned real roots (that is, the ones that admit isolation). 
Namely, one computes some complex discs, each covering a single real root and no 
other roots of the polynomial p(x ). Then the isolated roots are approximated fast 
by means of some specialized root refinement algorithms. The record and nearly 
optimal Boolean complexity at this stage has been obtained in the algorithm of 

[Em], Em]. 

Presently we achieve progress, so far missing, at the former stage of the 
real root isolation. Our algorithm performs this stage and consequently approxi¬ 
mates all the simple and well conditioned real roots within the same asymptotic 
Boolean complexity bounds of the papers [PT13] , [PT151 (see our Theorem [7]). 

We have also extended our algorithm to the isolation of complex, possi¬ 
bly multiple, roots and root clusters staying within the same (nearly optimal) 
asymptotic Boolean cost bound (see Section [6]). 

As in the papers PT13I . |PT15] , we approximate only simple and well con¬ 
ditioned real roots, but not the multiple and ill conditioned roots. The latter 
roots little affect the complexity of our algorithm, e.g., our cost estimate do not 
depend on the minimal distance between the roots and do not include the terms 
like log(Discr(p) -1 ). 

Our overall cost bound also matches the nearly optimal one of the papers 
[P95] and IP02] . Their algorithm approximates all complex roots of a polynomial, 
but combines a number of advanced techniques. This makes it much harder to 
implement and even to comprehend than our algorithm, which is much less 
involved, more transparent and more accessible for the implementation. 

We have tested our algorithm applied with the IEEE standard double pre¬ 
cision to some benchmark polynomials with small numbers of real roots. The 
test results are quite encouraging and are in good accordance with our formal 
study. The overall Boolean complexity of our isolation algorithm is dominated at 
the stage of performing the Dandelin’s auxiliary root-squaring iterations, whose 
Boolean cost grows fast when their number increases, but in our tests this num¬ 
ber grew very slowly as we increased the degree of the input polynomials from 
64 to 1024. 

Technically, we achieve our progress by incorporating the old algorithm of 
[S82] , for the approximation of the root radii, that is, the distances of the roots to 
the origin (we refer the reader to [0401 , G72) , |H74] , and [ B79] on the preceding 
works related to that algorithm). We feel that Schonhage in }S82] used only a 
small part of the potential power of the algorithm. Namely he applied it to the 










































rather modest task of the isolation of a single complex root, whereas we apply 
this simple but surprisingly efficient tool to the isolation of all simple and well 
conditioned real roots. We hope that this tool will be efficiently incorporated 
into other root-finders as well. 

Our another basic sub-algorithm performs multi-point evaluation of the poly¬ 
nomial p(x). The algorithm of (MB72I , recently used in the root-finders of 1PT13] , 
|PT14b) . |PT15j . and (SM15) . solves this problem at a low Boolean cost, but is 
numerically unstable and only works with extended precision. Performed with 
double precision, it produces corrupted output already for polynomials of mod¬ 
erate degree (say, about 50 or so). 

In some cases we can avoid this drawback, by applying the recent alternative 
numerically stable algorithms of IP15J and )Paj , whose Boolean cost matches the 
one of |MB72I as long as the outputs are required within the relative approxima¬ 
tion error bound l/2 b , for b = 0(log(n). This is certainly the case at the initial 
stage of our algorithm, when we compute the sign of p(x ) at 2 n points. 

We organize our presentation as follows. In the next section we cover some 
auxiliary results. In Section [3] we recall the Boolean complexity estimate for the 
approximation of the root radii. In Section[|]we describe our real root algorithm. 
In Section[5]we demonstrate it by working example. In Section[G]we extended our 
real root-finding algorithm to isolation of complex roots and root clusters. Sec¬ 
tion 0 the contribution of the second author, covers the results of our numerical 
tests. In Section [ 8 ] we briefly comment on further extension of our work. 

2 Some Definitions and Auxiliary Results 

Hereafter “flop” stands for “arithmetic operation”. “lg” stands for “log 2 ”. 

Ob(-) and Ob(-) denote the Boolean complexity up to some constant and 
poly-logarithmic factors, respectively. 

r is the overall bit-size of the coefficients of a polynomial p(x) of dTJ. 

Theorem 1. Count the roots of a polynomial with their multiplicity. Then a 
polynomial p(x) has an odd number of roots in the real line interval (a, /3) if and 
only if p(a)p((3) < 0 . 

Definition 1. D(z,p) = {a: : \x— z\ < p} denotes the closed disc with a complex 
center z and a radius p. Such a disc is 7 -isolated, for 7 > 1, if the disc D(z , 7 /a) 
contains no other roots of the polynomial p(x) of equation f71) . A root Xj is 7 - 
isolated if the disc D(xj, (7 — l)|^j|) contains no other roots of the polynomial 
p{x) besides Xj. 

The following theorem states that Newton’s iterations 

y 0 = z, yW = yW - p{yW)/pf(y™), h = 0,1 ,... ( 2 ) 

converge with quadratic rate globally, that is, right from the start, if they are ini¬ 
tialized in a 3 (n— l)-isolated disc containing a single simple root of a polynomial 
p(x) of 0 . 
















Theorem 2. Assume Newton’s iteration P|) for a polynomial p = p{x) of PJ1 
and let 0 < 3(n — l)|j/o — *i| < |yo — %j | for j = 2,... ,n. Then |yk — x\\ < 
2\yo - ^i|/2 2, ° for k = 0,1 ,.... 


Proof. This is T98 l Theorem 2.4], which strengthens [R87l Corollary 4.5]. 

The following two theorems state some upper bounds on the Boolean cost of 
certain fundamental polynomial computations. 

Theorem 3. Multi-point Polynomial Evaluation. Assume a real b > 1, a poly¬ 
nomial p(x) of pp, and k complex points Zi,...,Zk such that k > n. Write 
l = lg(max ™ =1 \zj\). Then, at the Boolean cost < 3 b ((& + r)k + Ikn), one can 
compute the values Vj such that \Vj — p(zj)\ < l/2 b for j = 1 ,,k. 

This is [K98j Theorem 3.9], also proved in [vdH08| . fKSal . and |PT14a] (cf. 
I’ll la Lemma 21], All proofs boil down to estimating the Boolean cost of the 
algorithm of ]MB72j , which relies on recursive polynomial division techniques of 
(F72I and fast FFT-based polynomial division of [S72] . 

Remark 1. Clearly, the same asymptotic bound applies to the Boolean complex¬ 
ity of a single Newton’s iteration performed concurrently at m points. 

Theorem 4. Given a real b > 1, a polynomial p(x) of pp, and the complex 
centers Zj and the radii pj ofm discs Dj = D(zj, pj), j = 1 ,..., m, each covering 
a single simple root of the polynomial p(x), write L = max 1 JL 1 lg(\zj\ + pj). Then, 
at the Boolean cost O s((& + L)n + rn 2 ), one can compute m new inclusion discs 
for the same roots, with the radii decreased by factors of at least 2 b provided that 

(i) all the m roots and all the m centers z i,..., z m are real or 

(ii) all the m discs D \,..., D m are 7 -isolated for a constant 7 > 1. 

Part (i) has been proved in [PT13] based on the algorithms of [PL99] and 
|PMRQTh7] (cf. also |PT15] and [SM 1 5j ) . Part (ii) has been proved in |PT14b] 
based on the approximation of the power sums of the roots of a polynomial. 

3 Root Radii and Their Estimation 

Definition 2. List the absolute values of the roots of p(x) in the non-increasing 
order, denote them rj = \xj\ for j = 1 ,..., n, r\ > r^ > • • • > r n , and call them 
the root radii of the polynomial p{x). 

The following result bounds the largest root radius 77 . 

Theorem 5. (See \VdS70y .) For a polynomial p{x) of m and ri = max " =1 \xj\, 
it holds that 

OArf /n< r\ < r(( for rf = 2 max \p n -ilPn\- 

i= 1 


( 3 ) 




































Remark 2. The theorem bounds all the root radii, but the paper IS05I (extend¬ 
ing the previous work in L798I . |K861 . and }H98j f presents stronger upper esti¬ 
mates for the positive roots of a polynomial p(x). By applying these estimates 
to the reverse polynomial p rev (x) = x n p(l/x) and to the polynomials p(—x) and 
p lev (—x), one can bound the positive roots from below and the negative roots 
from both below and above. 

Theorem 6 . Assume a polynomial p = p(x) of (if) and a positive A. Then, 
within the Boolean cost bound Ob{tti 2 ), one can compute approximations fj to 
all root radii rj such that 1/A < fj/rj < A for j = 1 provided that 

lg( ) = 0 (lg(n)), that is, \rj/rj — 1 | < c/n d for a fixed pair of constants 
c > 0 and d. 

This is [S82; Corollary 14.3]. At first the root radii are approximated at a 
dominated cost in the case of A = 2 n. See some details of the algorithm in 
[POOl Section 4]. [S82] Section 14] and [POOl Section 4] cite the related works 
10401 . 1072] , IH741 . 15851 . and one can also compare the relevant techniques of 
the power geometry in [B79] and 1B98I , developed for the study of algebraic and 
differential equations. 

In order to extend Theorem [ 6 ] to the case of A = (2ro ) 1 / 2 for any positive 
integer k , at first apply k Dandelin’s root-squaring iterations to the monic poly¬ 
nomial qo(x) = p(x)/p n (cf. m), that is, compute recursively the polynomials 

n 

qi{x) = {-l) n qi-i(y/x )qi-i{~Vx ) = - x 2 ), for i = 1,2,... (4) 

j=i 

Then approximate the root radii of qk(x) by applying Theorem[G]for A = 2 n 
and for p(x) replaced by qk{x). Finally approximate the root radii rj of p(x) as 
Vj = (rj fc '*) 1 / 2 . This enables us to decrease A from 2 n to at most 1 + c/n d = 
1 + 2 0 ( lg ( ra )) for any fixed pair of constants c > 0 and d by using k = 0 (lg(n)) 
Dandelin’s iterations. The cost bound of Theorem [G] follows. 

Remark 3. By using k = |"ln(s lg(n))] Dandelin’s root-squaring iterations, for 
s > 1 , the same algorithm supporting Theorem [ 6 ] approximates the root radii 
within the relative error bound - ln ^ 2rt ) ■ 

Corollary 1. Assume a real b > 1, a polynomial p(x) of (if, and a complex z. 
Then, within the Boolean cost bound Ob((t + n( 1 + /3))n 2 ), for (3 = lg(2 + |;s|), 
one can compute approximations fj ~ fj to the distances fj = \z — Xj | from the 
point z to all roots Xj of the polynomial p(x) such that 1/A < fj/fj < A, for 
j = 1 ,... ,n, provided that Ig(wzr) = 0 (lg(n)). 

Proof. Given a polynomial p(x) of (JT]) and a complex scalar z, one can compute 
the coefficients of the polynomial q(x) = p(x + z) by using 0(nlg(n)) flops (cf. 
[P01] b The root radii of the polynomial q{x) for a complex scalar z are equal to 
the distances \xj — z\ from the point z to the roots Xj of p(x). 


































Now let r q denote the bit-size of the coefficients of q(x), and let fj for j = 
1 ,... ,n denote its root radii, listed in the non-increasing order. Then, clearly, 
fj < r-j + \z\ for j = 1 , ...,n. Furthermore represent the coefficients of the 
polynomial q{x) with the bit-size r+n by using representation of p{x) with bit- 
size 0(T + n(l + /3)) for (3 = lg(2+|z|). By applying Theorem[6]to the polynomial 
q(x), extend the cost bounds from the root radii to the distances. 

4 Isolation and Approximation of Simple and 
Well-Conditioned Real Roots 

Algorithm 1. Approximation of simple and well-conditioned real roots. 

Input: two positive integers n and r such that 0 < r < n, a positive tolerance 
bound t, three real constants b, c and d such that b > 1 and c > 0 , the coefficients 
of a polynomial p{x) of equation (Q]j, and an upper estimate r' for the unknown 
number of its real roots, r' < n. 

Output: Approximations within the relative error bound l/2 b to some or all 
real roots Xi ,..., x r+ of the polynomial p(x), including all real roots that are 
both single and (1 + c/n d )-isolated. 

Computations: 

1. Compute approximations fi,...,f n to the root radii of a polynomial p(x) 
of W within the relative error bound c/n d , f\ < ••• < f n (see Theorem 
03 ). (Each approximation fi defines at most two real line intervals that can 
include real roots lying near if*. Overall we obtain 2 n candidate inclusion 
intervals T \,..., I 2 n, some of which can overlap or coincide with each other.) 

2. Define T = {ti,...,t n >}, the set of the distinct endpoints of these intervals 
in non-decreasing order, for n' < 4n. Compute the sign of the polynomial 
p{x) on this set. If the sign changes at two consecutive real points, t' = tj 
and t" = tj. (_i, select the line interval (t',t"). 

3. Apply the algorithm of IPTISf . \PT15( . which supports part (i) of Theorem ^) J 
concurrently to all selected line intervals. As soon as the algorithm decreases 
the length of a selected interval to or below l/2 b—1 , output the midpoint as 
an approximation to a real root. If r' such approximations have been output 
or if the cost of performing the iterations of the algorithm of IPTISf . \PI15A 
reaches the tolerance bound t, then stop the computations. (At stopping, there 
can remain some selected intervals of length exceeding l/2 h . We can output 
them with the label “intervals with real roots, apparently ill-conditioned”.) 

Compare the error bound on the root radii at Stage 1 with the definition of 
a single (1 + c/n d )-isolated real root of p(x) in Definition [T] and conclude that 
one of the inclusion intervals Tj, for 1 < j < 2 n, contains this root and no other 
roots of p{x). Now the correctness of the algorithm follows from Theorem [L] 
The Boolean cost at Stage 1 is 0 _b((t + l)n 2 ), for l = max^j lg(|xj|) (by 
virtue of Theorem [G]), at Stage 2 is Os(rn + In 2 ) (by virtue of Theorem 03 for 










6 = 2), and at Stage 3 is Ob((t + l)n 2 + bn) (by virtue of part (i) of Theorem 
2] for L = 0(1)). In view of ([3]), we can write l = 0 (t), and then the overall 
asymptotic cost bound, dominated at Stage 3, turns into Ob(tu 2 + bn). This is 
the same cost bound as at Stage 1 if 6 = 0(rn) and is the same as Stage 2 if in 
addition r = 0(1). Summarizing we obtain the following estimate. 

Theorem 7. Suppose that we are given the coefficients of a polynomial p(x) of 
equation m and three real constants b, c, and d such that b > 1 and c > 0. 
Then, at the cost Ob(tu 2 + bn), we can approximate all the (1 + c/n d )-isolated 
real roots of p(x) within the relative error bound 1 / 2 6 . 

Remark f. Before we begin the computations of Stage 2 of of Algorithm Q] we 
can narrow the ranges for positive and negative roots of the polynomial p(x) by 
following the recipes of Remark [2] Then at Stage 2 we can avoid the evaluation 
of the polynomial at the points lying outside these ranges. 

Remark 5. Before we started Stage 3, we can recompute the root radii for a 
smaller value c/n d and then check if some intervals, selected under this value, he 
strictly inside the old ones, sharing no endpoints, and hence are 7 -isolated for 
some 7 > 1. To such intervals we can apply a simplified version of the algorithm 
of [PTT 3 ] . [Em], which combines bisection and Newton’s iteration and does 
not involve the more advanced algorithm of double exponential sieve (also called 
the bisection of the exponent). 

Remark 6. For multi-point evaluation of the polynomial p(x), we can apply the 
algorithm of [MB72] with extended precision, but if the output is only required 
with a low precision, then instead we can apply the algorithms of |P15| and 
|Pa| with double precision at the same asymptotic Boolean cost. In particular at 
Stage 2 of Algorithm [l] we only seek a single bit (of the sign of p(x)) per point 
of evaluation, and at Stage 3 at the initial Newton’s iterations, crude values 
of p(x) can be sufficient. At Stage 2 of contracting r suspect real intervals, we 
perform order of 0(r lg( 6 )) evaluations, and if r is small, which is quite typically 
the case in the applications to algebraic and geometric computations, then it 
can be non-costly even to apply the so called Horner’s algorithm. It uses 2 nr 
arithmetic operations, that is, less than the algorithms of |MB72j . [ PI5] , and 
(Pa) , if r = o(log 2 (n)). 

Remark 7. If we only need to approximate the very well isolated real roots, 
we can apply Algorithm [I] at a lower cost for a larger isolation ratio. If we do 
not know how well the roots are isolated, we can at first apply Algorithm [T] 
assuming a larger ratio and, if the algorithm fails, re-apply it recursively, eah 
time assuming stronger isolation. 

5 Working Example 

Consider the following polynomial 

p(x) = 8x 7 + 16a; 6 + 16a; 5 + 16a; 4 — 23a; 3 — 30a; 2 + 3a; + 4. 














This product of the 4th degree Clrebyslrev polynomial of the first kind and the 
cubic polynomial x 3 + 2a; 2 + 3a; + 4 has five real roots and two complex conjugate 
non-real roots. All seven roots are well-separated, namely, 


0.3827 + 0.0000*, -0.3827 + O.OOOOi, 0.9239 + 0.0000a, -0.9239 + 0.0000a, 
-0.1747 + 1.5469a, -0.1747- 1.5469a, -1.6506 + 0.0000a. 

It took 14 root-squaring iterations to estimate the seven root-radii with a 
precision of at least 0.001: 

[0.3827,0.3827], [0.3826,0.3828], [0.9237,0.9240], [0.9238,0.9241], 

[1.5565,1.5570], [1.5565,1.5570], [1.6504,1.6507], 

The range [1.5565,1.5570] was for two roots. Taking the negations into ac¬ 
count, we obtained 12 intervals containing all real roots, among which we counted 
the intervals [—1.5565,-1.5570] and [1.5565,1.5570] twice: 

[-1.6507, -1.6504], [-1.5570, -1.5565], [-0.9241, -0.9238], 

[-0.9240, -0.9237], [-0.3828, -0.3826], [-0.3827, -0.3827], [0.3827,0.3827], 
[0.3826,0.3828], [0.9237,0.9240], [0.9238,0.9241], 

[1.5565,1.5570], [1.6504,1.6507], 

Some of the ranges for the remaining root radii overlapped pairwise, but 
our algorithm evaluated the polynomial p[x) at the endpoints of all the twelve 
intervals. Wherever there was the change of the sign, we concluded that the 
interval contains a root of the polynomial p(x). We obtained nine (rather than 
five or even seven) intervals with the sign changes, namely: 

[-1.6507, -1.6504] [-0.9241, -0.9238] [-0.9240, -0.9237] [-0.3828, -0.3826] 
[-0.3827, -0.3827] [0.3827,0.3827] [0.3826,0.3828] [0.9237,0.9240] [0.9238,0.9241] 

Then our algorithm initiated Newton’s iteration at the mid-point of each of 
these intervals. Whenever the Newton’s map fell outside the interval, the algo¬ 
rithm switched to bisection method instead of Newton’s. Very soon the algorithm 
converged to all real roots of p(x). 

The first Newton’s iteration produced the following real root estimates: 

-1.65062921, -0.92387954, -0.92387954, -0.38268343, 
-0.38268343,0.38268343,0.38268343,0.92387954,0.92387954. 

The second iteration produced the estimates 
-1.65062919, -0.92387953, -0.92387953, -0.38268343, -0.38268343, 
0.38268343,0.38268343,0.92387953,0.92387953. 

The maximum difference of the estimates output by the second and the 
third iterations was less than 10 -8 . Thus the algorithm stopped and returned 
the values of five distinct real roots of the polynomial p(x): 

-1.65062919, -0.92387953, -0.38268343,0.38268343,0.92387953. 



6 Isolation of Simple and Well-Conditioned Complex 
Roots by Means of Root-radii Approximation 

Algorithm 2. Isolation of Simple and Well-Conditioned Complex Roots. 
Input: two positive scalars p and e and the coefficients of a polynomial p(x) of 

W- 

Output: Set of approximations to the roots of the polynomial p(x) within p\J 2 
such that with probability at least 1 — e this set includes all roots having no other 
roots of the polynomial p{x) in their A 1 -neighborhoods for A' = (22.63 n 4 + 
2e)p/e. 

Initialization: Fix a reasonably large scalar p, say, p = 100. Generate a ran¬ 
dom value <f> under the uniform probability distribution in the range [ 71 / 8 , 37 t/ 8 ] . 
Computations: 

1. (Three Shifts of the Variable.) Compute the value rf = 2max" =1 \p n -i/p n \ of 

W- Then compute the coefficients of the three polynomials q(x) = p(x—prf), 
q-(x) = p(x - pr[]y/~A), and q</,(x) = p{x - pr[] . 

2. Compute approximations to all the n root radii of each of these three poly¬ 
nomials within the relative error bound p/((r/~ + 1 )p). [This defines three 
families of large thin annuli. Each family consists of n annuli and each an¬ 
nulus contains a root ofp(x). Multiple roots define multiple annuli. Clusters 
of roots define clusters of overlapping annuli.] 

3. Compute the intersections of all pairs of the annuli from the first two families 
and of the disc Z?(0,r/~). [We only care about the roots of p(x), and all of 
them lie in the disc D(0,r][). The intersection of each annulus with the 
disc D(0,r][) is close to a rectangular (vertical or horizontal) band of width 
at most p on the complex plane because every annulus has width at most 
P "A {ti + l )/?7 where we have chosen p large enough. The intersection of 
any pair of annuli from the two families is close to a square, with nearly 
vertical and nearly horizontal edges of length at most p[. The disc D(0,r][) 
contains a grid made up ofn 2 such squares, said to be nodes. [Some of them 
can be merged together (in the case of multiple roots) or nearly merged (in 
the case of root clusters).] 

f. For each annulus of the third family, determine its intersection with the nodes 
of the grid, and if the annulus intersects only a single node, then output its 
center and stop. 

Let us prove correctness of the algorithm. At first readily verify the following 
lemma. 

Lemma 1. Suppose that a line passes through a disc D(z,p') in the direction 
chosen at random under the uniform probability distribution of its angle in the 
range [a,a + 7 ] for 0 < 7 /( 271 ) < 1 . Then the line intersects a disc D(z',p') with 
a probability at most P = ^ arctan( |_/^,| )■ 

Next apply the lemma to a pair of nodes of the grid with distance A = \z — z'\ 
between their two centers z and z'. In this case 7 = 1/8 and the two nodes lie 



in the two discs D{z,p') and D(z',p' ) for p' = p\J 2. Furthermore assume that 
A p' , so that p^TTj ~ arctan(p^Tj) < 1.001 |-^ z ,i • Then the lemma implies 
that 

P < 16-\/2arctan(///Z\) < 22.6275 p/A. (5) 

Theorem 8. Let the grid of Algorithm [H have N p nodes overall, N p < n 2 . 
Suppose that an annulus of the third family intersects a node whose center has no 
centers of the other nodes of the grid within the distance A = 22.63(N p — 1 )p/e, 
for a fixed positive tolerance e. Then this annulus intersects another node of the 
grid with a probability at most e. 

Proof. Apply bound ([5]) to all pairs of the nodes of the grid. 

Now correctness of the algorithm follows because every root of the polynomial 
p(x) lies in some annulus of each of the three families. 

The Boolean complexity bounds for Stage 1 of Algorithm |T] are immediately 
extended to cover the Boolean complexity of Stages 1 and 2 of Algorithm [2] as 
long as log((rj f + l)rj/p = 0(log(n)) (cf. Corollary [T| . The Boolean complexity 
at Stage 3 of Algorithm |T] is dominated if we check separately whether the real 
parts and the imaginary parts of nodes overlap pairwise. Likewise the Boolean 
complexity at Stage 3 of Algorithm |T] is dominated if, for every annuls of the third 
family, we apply the bisection process, which begins with the median knot of the 
grid and recursively discards the knots that he above and below the annulus. 

Remark 8. We can modify Stage 2 of Algorithm |T] by collapsing every chain of 
m pairwise overlapping or coinciding root radii intervals (where m < n) into 
a single interval of relative width at most mp/((r^ + 1 ) 77 ) and by assigning 
to this interval multiplicity m. Such extended root radii define thicker annuli, 
of multiplicity m > 1. The intersection of two annuli of two families, having 
multiplicity mi and mi, respectively, defines a node of the new grid, to which 
we assign multiplicity min{mi,m 2 }. We do not change Stage 4, except that an 
output node of multiplicity m contains exactly m roots of the polynomial p(x), 
each counted according to its multiplicity. 

Remark 9. Clearly the assumptions of part (i) of Theorem 0] are satisfied for the 
centers of the output nodes of Algorithm [l] and so we can extend the algorithm 
to refining the approximations to the well conditioned roots of the polynomial 
p(x) given by these centers. For every node output by the algorithm of the 
previous remark, we can compress its superscribing disc fast, according to part 
(ii) of Theorem 01 as long the assumptions of that part are satisfied. 

Remark 10. We have chosen absolutely large shift values in order to simplify our 
analysis, although in this case we must use higher computational precision. One 
can apply the same algorithm for smaller value of 77. This heuristic variation may 
work and may allow to use lower precision of computing. 

Remark 11. If the nodes were points and the annuli straight lines, one could 
compute the minimum nonzero angle (j> between all the leftmost and all the 





rightmost nodes of the grid lying above. Then by directing the third family 
of the annuli by the angle <j>/2, one could ensure deterministically that every 
annulus passes through exactly a single node of the grid. This recipe can be 
extended to the case of sufficiently thin annuli and small nodes. If the width of 
the annuli is too large, we can decrease it by using more root-squaring iterations, 
but it is not clear for which class of inputs we can still keep the claimed overall 
cost bound under this recipe. 

7 The Tests for Real Root-finding with Algorithm |Tj 

We tested a variant of our isolation Algorithm |T] described in Remarks [5] [7] 
We applied it to polynomials of degrees n = 64,128,256,512,1024 having no 
clustered roots. In this variant we simplified Stage 3: instead of the algorithm of 
[PT131 , [PTlb we applied three Newton’s iterations, initialized at the center of 
the isolation interval output by Algorithm [1] In the case of poor convergence, 
we performed a few bisections of the output interval of Algorithm [Tj and then 
applied three Newton’s iterations again. We performed all computations with 
double precision and estimated the output errors by comparing our results with 
the outputs of MATLAB function ’’rootsQ”. 

We generated our input polynomials as the products of the Chebyshev poly¬ 
nomials of the first kind (having degree r and thus having r simple real roots) 
with polynomials of degree n-rof the following three families, 

1. Polynomials with real standard Gaussian random coefficients. 

2. Polynomials with complex standard Gaussian random coefficients. 

3. Polynomials with integral consecutive coefficients starting at 1, i.e., p(x) = 
1 + 2x + ... + (n — l)x n . 

Table [l] displays the number of Dandelin’s root-squaring iterations (in the col¬ 
umn “Iter”) and the maximum error (in the column “Error”). The tests showed 
very slow growth of the number of iterations, and consequently of the computa¬ 
tional cost of Stage 1, as n increased from 64 to 1024. 

8 Conclusions 

In our future study we plan to elaborate upon the extension of our isolation of 
multiple roots and root clusters, by incorporating the techniques that extend 
part (ii) of Theorem [I] and complementing them by the algorithm of |K98] . 
These techniques should enable us to split out the factor of p(x) whose root set 
is precisely the cluster of the roots of the polynomial p{x) lying in the output 
node, and then we can work on root-finding for the remaining complementary 
factor. Developing this approach, we are going to explore potential advantages 
of shifting the origin into the center of gravity of the roots, —p„_i/(np n ) (cf. 
[S82] ). or into the roots of higher order derivatives of p(x) (cf. [ MP13 , Section 

!5])- 

We plan to test our current and extended algorithms against the available 
real root-finders and to explore the chances for enhancing the efficiency of our 












Table 1. Real Root-finding with Using Root Radii Estimation 



Type 1 

Type 2 

Type 3 

E9 

m 

Uffl 

Error 

Q 

Error 

mu 

Error 
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B 

6.42E-11 

B 

6.42E-11 

4 

6.42E-11 



B 

3.16E-05 

B 

3.16E-05 

B 

3.16E-05 


m 


4.36E-03 


4.36E-03 

B 

4.36E-03 


m 

B 

6.42E-11 

11 

6.42E-11 

4 

6.42E-11 


m 

B 

3.16E-05 

B 

3.16E-05 

B 

3.16E-05 


m 


4.36E-03 

B 

4.36E-03 

B 

4.36E-03 


u 

B 

6.42E-11 

4 

6.42E-11 

4 

6.42E-11 


m 

B 

3.16E-05 

B 

3.16E-05 

B 

3.16E-05 


m 


4.36E-03 


4.36E-03 

B 

4.36E-03 


u 

B 

6.42E-11 

B 

6.42E-11 

4 

6.42E-11 



B 

3.16E-05 

B 

3.16E-05 

B 

3.16E-05 


EB 


4.36E-03 

B 

4.36E-03 

B 

4.36E-03 


H 

B 

6.42E-11 

B 

6.42E-11 

4 

6.42E-11 



B 

3.16E-05 

B 

3.16E-05 

B 

3.16E-05 


m 

B 

4.36E-03 

B 

4.36E-03 

B 

4.36E-03 


algorithms by means of their combination with some of the other highly efficient 
techniques known for root-finding. 
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